1 Introducción

En este trabajo se analiza la estructura multicapa de la red de transporte de la ciudad de Madrid a partir de dos conjuntos de datos reales proporcionados:

  • Un fichero con las paradas y estaciones de distintos modos de transporte.
  • Un fichero con información de viajes mensuales por distrito.

Los objetivos concretos son:

  • Construir una red multicapa donde cada capa representa un modo de transporte (bus, metro, bicimad, parking).
  • Estudiar la centralidad de las zonas de intercambio en cada capa.
  • Medir la versatilidad de las zonas (cómo cambia su importancia entre modos).
  • Realizar una visualización multicapa en la que se vean explícitamente las distintas capas conectadas.
  • Explorar de forma sencilla la vulnerabilidad del sistema ante fallos en una de las capas.

Toda la parte empírica se basa exclusivamente en los ficheros:

  • transporte_madrid_consolidado.csv
  • viajes_madrid.csv

No se utilizan ejemplos externos (como StarWars) para que el análisis esté completamente centrado en los datos de Madrid.

2 Carga y exploración de datos

library(tidyverse)
library(janitor)
library(sf)
library(dbscan)
library(FNN)
library(igraph)
library(scales)
library(knitr)

2.1 Lectura de los ficheros

Se asume que los ficheros CSV están en el mismo directorio que este documento RMarkdown. Si no es así, ajustar las rutas.

ruta_transporte <- "transporte_madrid_consolidado.csv"
ruta_viajes     <- "viajes_madrid.csv"

transporte_raw <- readr::read_delim(
  ruta_transporte,
  delim = ";",
  col_types = cols()
)

viajes_raw <- readr::read_delim(
  ruta_viajes,
  delim = ";",
  col_types = cols()
)

glimpse(transporte_raw)
## Rows: 14,188
## Columns: 6
## $ stop_name         <chr> "                                                   …
## $ stop_lat          <dbl> 40.37835, 40.38341, 40.36813, 40.37773, 40.40753, 40…
## $ stop_lon          <dbl> -3.75458, -3.62481, -3.62283, -3.67311, -3.59176, -3…
## $ transport_mode    <chr> "bicimad                                            …
## $ stop_id           <chr> "382 -Parque Eugenia de Montijo, 2                  …
## $ distrito_asignado <dbl> 2807911, 2807918, 2807918, 2807913, 2807919, 2807911…
glimpse(viajes_raw)
## Rows: 8,064
## Columns: 11
## $ fecha             <chr> "2024-01", "2024-01", "2024-01", "2024-01", "2024-01…
## $ name              <chr> "Centro                                            "…
## $ edad              <chr> "0-25", "0-25", "0-25", "0-25", "0-25", "0-25", "0-2…
## $ sexo              <chr> "hombre", "hombre", "hombre", "hombre", "mujer", "mu…
## $ numero_viajes     <chr> "0", "1", "2", "2+", "0", "1", "2", "2+", "0", "1", …
## $ personas          <dbl> 90541756, 18809207, 73248991, 187551534, 75357987, 1…
## $ Codigo            <dbl> 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, …
## $ Provincia         <chr> "Madrid                ", "Madrid                ", …
## $ distrito          <chr> "2807901 ", "2807901 ", "2807901 ", "2807901 ", "280…
## $ zona_pernoctacion <dbl> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ poblacion         <dbl> 140644, 140644, 140644, 140644, 140644, 140644, 1406…

2.2 Descripción básica de transporte_madrid_consolidado

transporte <- transporte_raw %>% 
  clean_names() %>% 
  mutate(
    transport_mode = str_trim(transport_mode),
    transport_mode = str_extract(transport_mode, "\\w+")
  )

table(transporte$transport_mode)
## 
## bicimad     bus   metro parking 
##     626   12430    1050      82
summary(select(transporte, stop_lat, stop_lon, distrito_asignado))
##     stop_lat        stop_lon      distrito_asignado
##  Min.   :40.28   Min.   :-3.875   Min.   :2807901  
##  1st Qu.:40.40   1st Qu.:-3.711   1st Qu.:2807905  
##  Median :40.43   Median :-3.689   Median :2807909  
##  Mean   :40.43   Mean   :-3.683   Mean   :2807910  
##  3rd Qu.:40.46   3rd Qu.:-3.656   3rd Qu.:2807915  
##  Max.   :40.56   Max.   :-3.447   Max.   :2807921

Variables clave:

  • stop_name: nombre de la parada/estación.
  • stop_lat, stop_lon: coordenadas geográficas.
  • transport_mode: modo de transporte (bus, metro, bicimad, parking).
  • stop_id: identificador de parada.
  • distrito_asignado: código de distrito (INE).

Esto define los nodos base del sistema de transporte.

2.3 Descripción básica de viajes_madrid

library(readr)
library(dplyr)
library(janitor)
library(stringr)

# Ajusta la ruta a tu fichero si hace falta
viajes <- read_delim("viajes_madrid.csv",
                     delim = ";",
                     locale = locale(decimal_mark = ","),
                     show_col_types = FALSE) %>% 
  clean_names()

# Comprobamos que personas es numérica y razonable
summary(viajes$personas)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4515   37486  110989  140469  186820  737342

Variables clave:

  • fecha: mes de referencia.
  • name: nombre del distrito.
  • edad, sexo: características sociodemográficas.
  • numero_viajes: categoría de nº de viajes (0, 1, 2, 2+).
  • personas: número de personas en esa celda.
  • distrito: código de distrito.
  • poblacion: población total del distrito.

2.3.1 Indicador de intensidad de viajes por distrito

Construimos un indicador sencillo de viajes medios por persona en cada distrito. Primero pasamos la categoría numero_viajes a un valor numérico aproximado y luego calculamos una media ponderada por el número de personas.

2.3.2 Indicador de viajes medios diarios por persona y distrito

La variable personas se encontraba codificada con formato decimal europeo, por lo que fue necesario convertirla explícitamente a formato numérico. Asimismo, la variable numero_viajes es categórica (0, 1, 2, 2+), por lo que se transformó en una aproximación numérica asignando el valor 3 a la categoría “2+”. El indicador final de viajes medios por persona se calculó como una media ponderada por el número de personas en cada categoría, obteniéndose valores en torno a 2 viajes diarios por persona, coherentes con patrones conocidos de movilidad urbana.

library(dplyr)
library(knitr)

viajes_distrito <- viajes %>% 
  mutate(
    # Pasamos las categorías de nº de viajes a un valor numérico aproximado
    viajes_cat = dplyr::case_when(
      numero_viajes == "0"  ~ 0,
      numero_viajes == "1"  ~ 1,
      numero_viajes == "2"  ~ 2,
      numero_viajes == "2+" ~ 3,
      TRUE ~ NA_real_
    )
  ) %>% 
  group_by(distrito) %>% 
  summarise(
    # Total de personas (peso) en el distrito
    personas_totales = sum(personas, na.rm = TRUE),
    # Total de viajes ponderados
    viajes_totales   = sum(viajes_cat * personas, na.rm = TRUE),
    # Viajes medios por persona
    viajes_medios_por_persona = viajes_totales / personas_totales,
    .groups = "drop"
  )

kable(
  head(viajes_distrito),
  digits = 3,
  caption = "Viajes medios diarios por persona y distrito"
)
Viajes medios diarios por persona y distrito
distrito personas_totales viajes_totales viajes_medios_por_persona
2807901 51888948 110673624 2.133
2807902 51542198 108276834 2.101
2807903 39251902 78164164 1.991
2807904 49039609 98112207 2.001
2807905 48349024 95741527 1.980
2807906 53679153 110757611 2.063
summary(viajes_distrito$viajes_medios_por_persona)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

El indicador de viajes medios diarios por persona presenta valores muy homogéneos entre distritos, situándose en un rango estrecho comprendido aproximadamente entre 1.98 y 2.13 viajes diarios por persona. Esta baja dispersión sugiere que, a escala distrital, los patrones de movilidad cotidiana en Madrid son relativamente uniformes, independientemente del volumen total de población o del número absoluto de desplazamientos. En particular, los distritos con mayor intensidad de viajes no necesariamente coinciden con los más poblados, lo que indica que la demanda de movilidad responde más a la estructura funcional del territorio y a la conectividad del sistema de transporte que al tamaño demográfico. Este resultado refuerza la idea de que la versatilidad multimodal y el papel estructural de ciertas zonas en la red pueden desempeñar un papel más relevante que la mera intensidad de viajes a la hora de caracterizar nodos clave del sistema de transporte urbano.

3 Construcción de la red multicapa

La idea es pasar de paradas individuales a zonas de intercambio (clusters de paradas cercanas), y luego construir una capa por modo de transporte.

3.1 1. Definición de zonas de intercambio (clustering espacial)

# Creamos objeto sf con coordenadas
transporte_sf <- transporte %>% 
  st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326, remove = FALSE) %>% 
  st_transform(25830)  # sistema métrico (UTM zona 30)

coords_m <- transporte_sf %>% 
  st_coordinates() %>% 
  as.matrix()

# Clustering DBSCAN: paradas a menos de 150 m se agrupan
set.seed(123)
db <- dbscan(coords_m, eps = 150, minPts = 2)

transporte_sf$cluster_id <- db$cluster

# Tratamos el ruido (cluster 0) como clusters individuales
transporte_sf <- transporte_sf %>% 
  mutate(
    cluster_id = if_else(
      cluster_id == 0,
      max(cluster_id) + row_number(),
      cluster_id
    )
  )

# Zona de intercambio = centroide del cluster
zonas <- transporte_sf %>% 
  st_drop_geometry() %>% 
  group_by(cluster_id) %>% 
  summarise(
    stop_lat = mean(stop_lat),
    stop_lon = mean(stop_lon),
    # distrito asignado más frecuente en el cluster
    distrito_asignado = as.integer(names(which.max(table(distrito_asignado)))),
    modos_disponibles = paste(sort(unique(transport_mode)), collapse = ","),
    .groups = "drop"
  )

nrow(zonas)
## [1] 1558
head(zonas)

Cada cluster_id representa una zona de intercambio multimodal.

Representamos la distribución espacial aproximada (proyección geográfica simple):

zonas_sf <- zonas %>% 
  st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326)

plot(zonas_sf["distrito_asignado"], main = "Zonas de intercambio por distrito")

3.2 2. Nodos por modo de transporte

modos <- sort(unique(transporte_sf$transport_mode))

zonas_modo <- transporte_sf %>% 
  st_drop_geometry() %>% 
  select(cluster_id, transport_mode) %>% 
  distinct() %>% 
  mutate(presente = 1) %>% 
  pivot_wider(
    names_from = transport_mode,
    values_from = presente,
    values_fill = 0
  )

zonas_attrs <- zonas %>% 
  left_join(zonas_modo, by = "cluster_id") %>% 
  mutate(across(all_of(modos), ~replace_na(., 0)))

kable(head(zonas_attrs), caption = "Atributos básicos de las zonas por modo")
Atributos básicos de las zonas por modo
cluster_id stop_lat stop_lon distrito_asignado modos_disponibles bicimad bus metro parking
1 40.38368 -3.624114 2807918 bicimad,bus,metro 1 1 1 0
2 40.38096 -3.728060 2807911 bicimad,bus,metro 1 1 1 0
3 40.38626 -3.719795 2807912 bicimad,bus,metro 1 1 1 0
4 40.37886 -3.732537 2807911 bicimad,bus 1 1 0 0
5 40.40855 -3.672909 2807903 bicimad,bus 1 1 0 0
6 40.40380 -3.706547 2807901 bicimad,metro 1 0 1 0

3.3 3. Unimos la información de viajes por distrito

# Forzamos el mismo tipo en ambas columnas de distrito
zonas_attrs <- zonas_attrs %>% 
  mutate(distrito_asignado = as.integer(distrito_asignado))

viajes_distrito <- viajes_distrito %>% 
  mutate(distrito = as.integer(distrito))

# Join seguro
zonas_attrs <- zonas_attrs %>% 
  left_join(
    viajes_distrito,
    by = c("distrito_asignado" = "distrito")
  )

summary(zonas_attrs$viajes_medios_por_persona)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

En muchas zonas puede que no haya dato (si el distrito no aparece o hay problemas de asignación), pero para el análisis de red esto no es un problema; simplemente tendremos algunos NA.

3.4 4. Construcción de redes intra-capa

Para cada modo:

  • Consideramos solo las zonas donde ese modo está presente.
  • Unimos cada zona con sus k vecinos más cercanos (aquí k = 3) basándonos en la distancia euclídea de los centroides.
  • Asignamos un peso inversamente proporcional a la distancia.
# Centroides en coordenadas métricas
zonas_centroides <- zonas_sf %>% 
  st_transform(25830) %>% 
  st_coordinates() %>% 
  as_tibble() %>% 
  rename(x = X, y = Y) %>% 
  mutate(cluster_id = zonas$cluster_id)

build_edges_modo <- function(modo, k = 3) {
  zonas_modo_ids <- zonas_attrs %>% 
    filter(.data[[modo]] == 1) %>% 
    pull(cluster_id)

  if (length(zonas_modo_ids) <= 1) {
    return(tibble(modo = character(), from = integer(),
                  to = integer(), dist = numeric(), w = numeric()))
  }

  coords_modo <- zonas_centroides %>% 
    filter(cluster_id %in% zonas_modo_ids) %>% 
    arrange(cluster_id)

  mat <- as.matrix(coords_modo[, c("x", "y")])
  k_eff <- min(k, nrow(mat) - 1)
  knn <- FNN::get.knn(mat, k = k_eff)

  from <- rep(coords_modo$cluster_id, each = k_eff)
  to   <- coords_modo$cluster_id[c(knn$nn.index)]
  dist <- as.numeric(knn$nn.dist)

  tibble(
    modo  = modo,
    from  = from,
    to    = to,
    dist  = dist,
    w     = 1 / (1 + dist)
  ) %>% 
    distinct() %>% 
    filter(from != to)
}

edges_intra <- map_df(modos, build_edges_modo)

kable(head(edges_intra), caption = "Ejemplo de aristas intra-capa")
Ejemplo de aristas intra-capa
modo from to dist w
bicimad 1 288 649.2619 0.0015378
bicimad 1 300 321.5842 0.0031000
bicimad 1 12582 476.5779 0.0020939
bicimad 2 283 506.9860 0.0019686
bicimad 2 927 343.8273 0.0029000
bicimad 3 29 398.1358 0.0025054

3.5 5. Grafos por capa (lista g_list)

nodos_ids <- zonas_attrs$cluster_id

build_graph_modo <- function(modo) {
  e_modo <- edges_intra %>% 
    filter(modo == !!modo) %>% 
    select(from, to, w)

  g <- graph_from_data_frame(
    d = e_modo,
    directed = FALSE,
    vertices = zonas_attrs %>% 
      mutate(name = cluster_id)
  )

  missing_nodes <- setdiff(nodos_ids, as.integer(V(g)$name))
  if (length(missing_nodes) > 0) {
    g <- g + vertices(as.character(missing_nodes))
  }

  g
}

g_list <- map(modos, build_graph_modo)
names(g_list) <- modos

g_list
## $bicimad
## IGRAPH 6bce25f UN-- 1558 1401 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bce25f (vertex names):
##  [1]  1--  288  1--  300  1--12582  2--  283  2--  927  3--   29  3--  735
##  [8]  3--  248  4--  372  4-- 6182  4--15480  5--  583  5--  778  5--  626
## [15]  6--  296  6-- 1121  6--  813  7--  239  7--  327  7--  370  8--  209
## [22]  8--  152  8--  145 11-- 1182 11--  813 11-- 6031 12-- 6031 12--  115
## [29] 12-- 4280 13--   47 13--  174 13-- 9493 14-- 1256 14--  621 14--  814
## + ... omitted several edges
## 
## $bus
## IGRAPH 6bd1215 UN-- 1558 4036 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1215 (vertex names):
##  [1]  1--1289  1--1026  1--1303  2-- 652  2-- 283  2--  29  3-- 819  3-- 390
##  [9]  3-- 844  4-- 767  4--1288  4-- 638  5-- 158  5-- 422  5-- 969  7--  20
## [17]  7--  18  7--  18  8-- 997  8-- 785  8-- 464  9-- 734  9-- 365  9-- 588
## [25] 10--1121 10--1011 10-- 239 11--  31 11--1118 11--  33 12--  32 12-- 365
## [33] 12-- 532 13-- 296 13-- 938 13--1255 14-- 854 14-- 599 14--1182 15--  43
## + ... omitted several edges
## 
## $metro
## IGRAPH 6bd1ffc UN-- 1558 682 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1ffc (vertex names):
##  [1]  1--230  1--229  1--  2  2--239  2--239  2--248  3--372  3--112  3--187
## [10]  6--187  6--248  6--917  7--287  7-- 68  7--145 11--753 11--220 11-- 49
## [19] 12--235 12-- 47 12--174 13--249 13-- 82 13--423 17--814 17-- 59 17-- 57
## [28] 18-- 57 18--516 18--564 24-- 30 24--446 24--223 27--222 27--437 27--118
## [37] 28--348 28--523 28--243 30--170 30--145 30--573
## + ... omitted several edges
## 
## $parking
## IGRAPH 6bd2c18 UN-- 1558 161 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd2c18 (vertex names):
##  [1]   7--   29   7--14152   7--  766   8--11008   8--  239   8--  245
##  [7]  27-- 1052  27--  814  27--  814  28--  446  28--  657  28--  164
## [13]  29--  573  29--  766  29--  941  41--  423  41--11008  41--  562
## [19]  29--   55  55--  423  55--  241  56--  939  56--  534  56--  678
## [25]  67--  613  67--14152  67--  630  71--  630  71--  814  71--  241
## + ... omitted several edges

4 Visualización multicapa con muxViz

En esta sección generamos la visualización clave: un gráfico 3D en el que se ven las capas (modos) apiladas, con los nodos alineados entre capas.

Para ello usamos la función plot_multiplex3D de muxViz. Se asume que muxViz está instalado; si no lo está, el código no fallará, pero mostrará un mensaje.

if (requireNamespace("muxViz", quietly = TRUE)) {
  library(muxViz)

  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  node_multimodalidad <- zonas_attrs %>% 
    mutate(
      modos_activos = rowSums(across(all_of(modos)))
    ) %>% 
    pull(modos_activos)

  node_size_values <- rescale(node_multimodalidad, to = c(0.5, 2))

  muxViz::plot_multiplex3D(
    g.list           = g_list,
    layer.colors     = layer_colors[names(g_list)],
    layer.labels     = "auto",     # <- clave: escalar, no vector
    node.size.values = node_size_values,
    node.size.scale  = 1.5,
    edge.size.scale  = 0.5,
    layout           = "fr",
    show.aggregate   = TRUE,
    show.nodeLabels  = FALSE
  )
} else {
  message("El paquete 'muxViz' no está instalado. Instálalo con devtools::install_github('manlius/muxViz') para ver la visualización 3D.")
}

Dado que la función plot_multiplex3D de muxViz presenta problemas en la instalación utilizada (bug interno con layer.labels), se ha implementado un gráfico multicapa alternativo usando tidygraph y ggraph, que representa las capas apiladas y las conexiones intra- e inter-capa de forma equivalente.

# Gráfico multicapa alternativo con tidygraph + ggraph

library(tidygraph)
library(ggraph)

# Orden de las capas y un índice para apilarlas verticalmente
layer_order <- modos
layer_index <- setNames(seq_along(layer_order), layer_order)

# Distancia vertical entre capas (en coordenadas proyectadas)
offset <- (max(zonas_centroides$y) - min(zonas_centroides$y)) * 1.2

# Nodos multilayer: mismo cluster_id repetido por capa, apilado en vertical
nodes_multi <- purrr::map_dfr(layer_order, function(m) {
  zonas_centroides %>%
    transmute(
      node       = paste(cluster_id, m, sep = "_"),  # id único por nodo-capa
      cluster_id = cluster_id,
      layer      = m,
      layer_idx  = layer_index[m],
      x          = x,
      y          = y + (layer_index[m] - 1) * offset
    )
})

# Aristas intra-capa: usamos edges_intra y las asociamos a cada capa
edges_intra_multi <- edges_intra %>%
  transmute(
    from  = paste(from, modo, sep = "_"),
    to    = paste(to,   modo, sep = "_"),
    type  = "intra",
    layer = modo
  )

# Aristas inter-capa: unen el MISMO cluster_id entre modos presentes (verticales)
inter_edges <- purrr::map_dfr(zonas_attrs$cluster_id, function(cid) {
  # modos presentes en esa zona
  fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
  modos_presentes <- modos[as.logical(unlist(fila))]
  
  if (length(modos_presentes) < 2) return(tibble())
  
  comb <- t(combn(modos_presentes, 2))
  tibble(
    from  = paste(cid, comb[, 1], sep = "_"),
    to    = paste(cid, comb[, 2], sep = "_"),
    type  = "inter",
    layer = "inter"
  )
})

edges_all <- bind_rows(edges_intra_multi, inter_edges)

# Grafo multilayer
g_multi <- tbl_graph(nodes = nodes_multi, edges = edges_all, directed = FALSE)

# Colores de capa (igual que antes)
layer_colors <- c(
  bicimad = "#1a9641",
  bus     = "#d7191c",
  metro   = "#2c7bb6",
  parking = "#fdae61"
)

ggraph(g_multi, layout = "manual", x = x, y = y) +
  geom_edge_link(aes(alpha = type, linetype = type), show.legend = TRUE) +
  geom_node_point(aes(color = layer), size = 1.2) +
  scale_edge_alpha_manual(values = c(intra = 0.3, inter = 0.8)) +
  scale_edge_linetype_manual(values = c(intra = "solid", inter = "dashed")) +
  scale_color_manual(values = layer_colors) +
  theme_void() +
  labs(
    title = "Red multicapa de transporte de Madrid",
    subtitle = "Capas (modos) apiladas verticalmente y conectadas por zonas de intercambio",
    color = "Capa",
    linetype = "Tipo de enlace",
    alpha = "Tipo de enlace"
  )

Este gráfico muestra claramente:

  • Las capas (bus, metro, bicimad, parking) como planos separados.
  • Los nodos (zonas de intercambio) alineados entre capas.
  • La estructura de conexiones en cada capa y la capa agregada.
# ===========================================
# Visualización 3D multicapa (inter-capa con scatter3d lines)
# ===========================================

if (requireNamespace("plotly", quietly = TRUE)) {

  library(plotly)
  library(dplyr)
  library(purrr)
  library(scales)
  library(tidyr)

  # ---- 1. Capas y eje Z
  layer_order <- modos
  layer_index <- setNames(seq_along(layer_order) - 1, layer_order)  # 0,1,2,3

  # ---- 2. Nodos multilayer: solo zonas donde ese modo existe
  nodes_multi_3d <- map_dfr(layer_order, function(m) {
    zonas_attrs %>%
      filter(.data[[m]] == 1) %>%                # solo zonas con ese modo
      select(cluster_id) %>%
      inner_join(zonas_centroides, by = "cluster_id") %>%
      transmute(
        node       = paste(cluster_id, m, sep = "_"),
        cluster_id = cluster_id,
        layer      = m,
        x_raw      = x,
        y_raw      = y
      )
  })

  # Normalizamos X, Y y asignamos Z
  nodes_multi_3d <- nodes_multi_3d %>%
    mutate(
      x = rescale(x_raw, to = c(0, 1)),
      y = rescale(y_raw, to = c(0, 1)),
      z = layer_index[layer]
    )

  # ---- 3. Multimodalidad por zona
  node_multimodalidad <- zonas_attrs %>%
    mutate(modos_activos = rowSums(across(all_of(modos)))) %>%
    select(cluster_id, modos_activos)

  # Clusters con al menos 2 modos
  clusters_interes <- node_multimodalidad %>%
    filter(modos_activos >= 2) %>%
    pull(cluster_id)

  # ---- 4. Aristas inter-capa: todos los pares de modos presentes
  edges_inter_3d <- map_dfr(clusters_interes, function(cid) {

    fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
    modos_presentes <- modos[as.logical(unlist(fila))]

    if (length(modos_presentes) < 2) return(tibble())

    comb <- t(combn(modos_presentes, 2))  # todos los pares posibles

    tibble(
      from       = paste(cid, comb[, 1], sep = "_"),
      to         = paste(cid, comb[, 2], sep = "_"),
      cluster_id = cid
    )
  })

  cat("Número de aristas inter-capa generadas:", nrow(edges_inter_3d), "\n")

  # ---- 5. Posiciones de nodos para las aristas
  pos_nodes <- nodes_multi_3d %>%
    select(node, x, y, z, layer)

  edges_plot <- edges_inter_3d %>%
    left_join(pos_nodes, by = c("from" = "node")) %>%
    rename(
      x  = x,
      y  = y,
      z  = z
    ) %>%
    left_join(pos_nodes, by = c("to" = "node"), suffix = c("", "_to")) %>%
    rename(
      xend = x_to,
      yend = y_to,
      zend = z_to
    ) %>%
    filter(
      !is.na(x), !is.na(y), !is.na(z),
      !is.na(xend), !is.na(yend), !is.na(zend)
    ) %>%
    select(x, y, z, xend, yend, zend, cluster_id)

  cat("Aristas con coordenadas válidas:", nrow(edges_plot), "\n")

  # ---- 6. Convertimos aristas a formato "long" para scatter3d lines
  # Para cada arista: (x,y,z) -> (xend,yend,zend) -> NA como separador
  if (nrow(edges_plot) > 0) {
    edges_long <- edges_plot %>%
      mutate(seg_id = row_number()) %>%
      mutate(
        x  = map2(x,  xend,  ~c(.x, .y, NA_real_)),
        y  = map2(y,  yend,  ~c(.x, .y, NA_real_)),
        z  = map2(z,  zend,  ~c(.x, .y, NA_real_))
      ) %>%
      select(seg_id, x, y, z) %>%
      unnest(c(x, y, z))
  } else {
    edges_long <- tibble(x = numeric(), y = numeric(), z = numeric())
  }

  cat("Puntos totales en edges_long:", nrow(edges_long), "\n")

  # ---- 7. Añadimos multimodalidad y tamaño de nodo
  nodes_multi_3d <- nodes_multi_3d %>%
    left_join(node_multimodalidad, by = "cluster_id") %>%
    mutate(
      modos_activos = replace_na(modos_activos, 1),
      size = rescale(modos_activos, to = c(4, 11))
    )

  # ---- 8. Colores por capa
  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  # ---- 9. Gráfico 3D
  p3d <- plot_ly()

  # Nodos
  p3d <- p3d %>%
    add_trace(
      data = nodes_multi_3d,
      type = "scatter3d",
      mode = "markers",
      x = ~x, y = ~y, z = ~z,
      color = ~layer,
      colors = layer_colors,
      marker = list(size = ~size),
      hoverinfo = "text",
      text = ~paste(
        "Zona:", cluster_id,
        "<br>Modo:", layer,
        "<br>Modos activos (total cluster):", modos_activos
      )
    )

  # Aristas inter-capa como líneas rojas gruesas
  if (nrow(edges_long) > 0) {
    p3d <- p3d %>%
      add_trace(
        data = edges_long,
        type = "scatter3d",
        mode = "lines",
        x = ~x, y = ~y, z = ~z,
        line = list(color = "black", width = 4),
        showlegend = FALSE,
        hoverinfo = "none"
      )
  } else {
    message("edges_long está vacío: no se han podido dibujar aristas inter-capa.")
  }

  p3d <- p3d %>%
    layout(
      title = "Red multicapa de transporte de Madrid (3D, enlaces entre capas)",
      scene = list(
        xaxis = list(title = "X (normalizada)"),
        yaxis = list(title = "Y (normalizada)"),
        zaxis = list(
          title = "Capa (modo de transporte)",
          tickvals = layer_index,
          ticktext = layer_order
        ),
        aspectmode = "cube"
      )
    )

  p3d

} else {
  message("Instala 'plotly' con install.packages('plotly') para ver la visualización 3D.")
}
## Número de aristas inter-capa generadas: 728 
## Aristas con coordenadas válidas: 728 
## Puntos totales en edges_long: 2184

Interpretación de la visualización 3D multicapa

La visualización tridimensional obtenida permite observar de forma explícita la estructura multicapa del sistema de transporte urbano de Madrid, donde cada plano horizontal representa un modo de transporte (bicimad, bus, metro y parking) y las conexiones verticales indican zonas de intercambio multimodal.

Separación clara de capas y estructura espacial

En primer lugar, se aprecia una separación nítida entre las cuatro capas, lo que confirma que cada modo de transporte presenta una distribución espacial propia:

  • La capa bicimad aparece más concentrada espacialmente, reflejando su carácter urbano y centralizado.
  • La capa bus muestra una cobertura más extensa, coherente con su función de red capilar que conecta barrios periféricos.
  • La capa metro se sitúa entre ambas, con una estructura más compacta que refleja la jerarquía de estaciones y líneas.
  • La capa parking aparece claramente diferenciada y menos densa, lo que concuerda con su función de infraestructura de apoyo más localizada.

Esta separación valida la decisión metodológica de modelar el sistema como una red multicapa, ya que una red agregada ocultaría estas diferencias estructurales.

Uniones verticales como indicador de multimodalidad

Las líneas verticales que conectan nodos entre capas representan zonas donde coexisten varios modos de transporte. Estas conexiones no aparecen de forma homogénea, sino que se concentran en un subconjunto reducido de zonas, lo que pone de manifiesto una fuerte heterogeneidad en la multimodalidad del sistema.

En términos de teoría de redes:

  • Estas zonas actúan como nodos altamente versátiles, ya que participan simultáneamente en varias capas.

  • Son puntos críticos para la integración funcional del sistema, permitiendo el transbordo eficiente entre modos.

  • Su posición estructural las convierte en puntos de control del flujo global de movilidad.

Visualmente, estas zonas destacan como “columnas” que atraviesan varias capas, reforzando la idea de que la conectividad global del sistema depende de un número limitado de nodos clave.

Versatilidad y roles multinodales

La figura ilustra de manera intuitiva el concepto de versatilidad de nodo:

  • La mayoría de los nodos aparecen en una sola capa, actuando como nodos especializados (por ejemplo, paradas de bus o estaciones de bici aisladas).
  • En contraste, un número reducido de nodos conecta tres o incluso cuatro capas, desempeñando un rol multinodal esencial.

Este patrón es consistente con la literatura sobre redes multicapas, donde se observa que la importancia de un nodo no puede evaluarse correctamente desde una sola capa. La visualización 3D confirma que los nodos versátiles no solo existen, sino que son estructuralmente visibles cuando se representan explícitamente las capas y sus interdependencias.

Implicaciones para vulnerabilidad y resiliencia

Desde el punto de vista de la vulnerabilidad del sistema, la visualización sugiere una conclusión relevante:

  • La conectividad multimodal del transporte madrileño depende de un conjunto relativamente pequeño de nodos altamente interconectados entre capas.
  • Una perturbación localizada en estos nodos (por ejemplo, cierre de un intercambiador o fallo simultáneo de varios modos en una zona) podría propagarse rápidamente al resto del sistema, afectando a múltiples capas.

Este comportamiento es coherente con los modelos de cascadas de fallo en redes interdependientes descritos en la literatura (Buldyrev et al., 2010; Duan et al., 2019), donde la interdependencia aumenta tanto la eficiencia como la fragilidad del sistema.

Valor añadido de la representación multicapa

Finalmente, esta visualización demuestra el valor añadido del enfoque multicapa frente a representaciones tradicionales:

  • Permite identificar de forma directa los intercambiadores críticos.
  • Hace visible la estructura de interdependencia entre modos.
  • Facilita la interpretación de conceptos abstractos como versatilidad, roles multinodales y vulnerabilidad sistémica.

En conjunto, el gráfico 3D no solo actúa como herramienta exploratoria, sino como una evidencia visual clara de que el sistema de transporte de Madrid funciona como una red interconectada de capas, cuya eficiencia y resiliencia dependen de la correcta gestión de sus nodos multimodales clave.

5 Centralidad y versatilidad de las zonas

5.1 1. Medidas de centralidad por capa

Calculamos grado y fuerza (suma de pesos) para cada zona en cada modo.

centralidad_por_capa <- map2_df(
  .x = g_list,
  .y = names(g_list),
  .f = function(g, modo) {
    tibble(
      cluster_id = as.integer(V(g)$name),
      modo       = modo,
      degree     = degree(g),
      strength   = strength(g, weights = E(g)$w)
    )
  }
)

kable(head(centralidad_por_capa), caption = "Centralidad (grado y fuerza) por capa")
Centralidad (grado y fuerza) por capa
cluster_id modo degree strength
1 bicimad 5 0.0096077
2 bicimad 4 0.0090672
3 bicimad 4 0.0078289
4 bicimad 6 0.0134323
5 bicimad 7 0.0140216
6 bicimad 6 0.0139194

Distribución del grado por modo:

centralidad_por_capa %>% 
  ggplot(aes(x = degree, fill = modo)) +
  geom_histogram(alpha = 0.6, bins = 30, position = "identity") +
  facet_wrap(~ modo, scales = "free_y") +
  theme_minimal() +
  labs(
    title = "Distribución del grado por modo de transporte",
    x = "Grado",
    y = "Frecuencia"
  )

5.2 2. Versatilidad de grado

Definimos una medida de versatilidad basada en la entropía de Shannon de la distribución de grado entre modos.

mat_grado <- centralidad_por_capa %>% 
  select(cluster_id, modo, degree) %>% 
  pivot_wider(
    names_from = modo,
    values_from = degree,
    values_fill = 0
  ) %>% 
  arrange(cluster_id)

grado_mat <- as.matrix(mat_grado[, modos])
rownames(grado_mat) <- mat_grado$cluster_id

calc_versatilidad <- function(vec) {
  if (sum(vec) == 0) return(NA_real_)
  p <- vec / sum(vec)
  p <- p[p > 0]
  H <- -sum(p * log(p))
  H / log(length(vec))
}

versatilidad <- apply(grado_mat, 1, calc_versatilidad)

versatilidad_df <- tibble(
  cluster_id = as.integer(names(versatilidad)),
  versatilidad_grado = as.numeric(versatilidad)
) %>% 
  left_join(zonas_attrs, by = "cluster_id")

summary(versatilidad_df$versatilidad_grado)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1477  0.0000  0.9983

5.2.1 Nodos más versátiles y más especializados

top_versatiles <- versatilidad_df %>% 
  arrange(desc(versatilidad_grado)) %>% 
  slice_head(n = 15) %>% 
  select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)

top_especializados <- versatilidad_df %>% 
  arrange(versatilidad_grado) %>% 
  slice_head(n = 15) %>% 
  select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)

kable(top_versatiles, caption = "Top 15 zonas más versátiles")
Top 15 zonas más versátiles
cluster_id versatilidad_grado modos_disponibles distrito_asignado
164 0.9983119 bicimad,bus,metro,parking 2807907
446 0.9976183 bicimad,bus,metro,parking 2807904
101 0.9968869 bicimad,bus,metro,parking 2807907
562 0.9963891 bicimad,bus,metro,parking 2807920
41 0.9955380 bicimad,bus,metro,parking 2807907
56 0.9949672 bicimad,bus,metro,parking 2807904
283 0.9949672 bicimad,bus,metro,parking 2807903
239 0.9946468 bicimad,bus,metro,parking 2807901
814 0.9934804 bicimad,bus,metro,parking 2807904
71 0.9927376 bicimad,bus,metro,parking 2807904
339 0.9927376 bicimad,bus,metro,parking 2807911
534 0.9927376 bicimad,bus,metro,parking 2807908
7 0.9926388 bicimad,bus,metro,parking 2807901
251 0.9892423 bicimad,bus,metro,parking 2807916
116 0.9863237 bicimad,bus,metro,parking 2807907
kable(top_especializados, caption = "Top 15 zonas más especializadas")
Top 15 zonas más especializadas
cluster_id versatilidad_grado modos_disponibles distrito_asignado
9 0 bus 2807903
10 0 bus 2807909
15 0 bus 2807921
17 0 metro 2807921
19 0 bus 2807921
20 0 bus 2807921
21 0 bus 2807921
22 0 bus 2807908
26 0 bus 2807911
31 0 bus 2807921
32 0 bus 2807921
33 0 bus 2807921
34 0 bus 2807911
36 0 bus 2807911
38 0 bus 2807920

Interpretación cualitativa (a redactar en la memoria):

  • Zonas con alta versatilidad: concentran paradas de varios modos y son relevantes en varias capas.
  • Zonas especializadas: la conectividad se concentra en un único modo.

5.3 3. Relación entre versatilidad e intensidad de viajes

versatilidad_df %>% 
  ggplot(aes(x = viajes_medios_por_persona, y = versatilidad_grado)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE) +
  theme_minimal() +
  labs(
    title = "Versatilidad de las zonas vs viajes medios por persona",
    x = "Viajes medios por persona (distrito)",
    y = "Versatilidad de grado"
  )

Comprobamos cuántas observaciones completas hay y calculamos la correlación solo si tiene sentido estadísticamente (al menos 3 observaciones completas):

datos_cor <- versatilidad_df %>% 
  filter(!is.na(versatilidad_grado),
         !is.na(viajes_medios_por_persona))

n_obs <- nrow(datos_cor)
n_obs
## [1] 1558
if (n_obs >= 3) {
  cor_test <- cor.test(
    datos_cor$versatilidad_grado,
    datos_cor$viajes_medios_por_persona
  )
  cor_test
} else {
  cat("No hay suficientes observaciones completas para estimar de forma fiable la correlación entre versatilidad y viajes.
")
}
## 
##  Pearson's product-moment correlation
## 
## data:  datos_cor$versatilidad_grado and datos_cor$viajes_medios_por_persona
## t = 0.13627, df = 1556, p-value = 0.8916
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04621554  0.05310757
## sample estimates:
##         cor 
## 0.003454533

Esto evita el error de “not enough finite observations” y hace el análisis más robusto.

6 Vulnerabilidad de la red de transporte

En lugar de implementar un modelo muy complejo de cascadas, realizamos un análisis sencillo pero informativo de robustez estructural:

  • Tomamos la red agregada (unión de todas las capas).
  • Eliminamos nodos de forma dirigida (los más centrales) en una capa y observamos cómo disminuye el tamaño de la componente gigante agregada.

6.1 1. Red agregada

# Matrices de adyacencia por capa
adj_list <- lapply(g_list, function(g) {
  as.matrix(as_adj(g, attr = "w", sparse = FALSE))
})

# Aseguramos mismo orden de nodos
orden_nodos <- order(as.integer(V(g_list[[1]])$name))
adj_list <- lapply(adj_list, function(m) m[orden_nodos, orden_nodos])

adj_agregada <- Reduce("+", adj_list)

g_agregado <- graph_from_adjacency_matrix(adj_agregada, mode = "undirected", weighted = TRUE)

comp <- components(g_agregado)
tamano_gcc_inicial <- max(comp$csize)
tamano_gcc_inicial
## [1] 1558

6.2 2. Robustez frente a fallos en una capa

Analizamos qué ocurre si vamos eliminando nodos con mayor grado en la capa metro y miramos el tamaño relativo de la componente gigante en la red agregada.

g_metro <- g_list[["metro"]]

deg_metro <- degree(g_metro)
orden_fallo <- names(sort(deg_metro, decreasing = TRUE))  # orden de fallo: más conectados primero

fracciones <- seq(0, 0.5, by = 0.05)

resultado_robustez <- map_df(fracciones, function(f) {
  n_remove <- ceiling(f * length(orden_fallo))
  nodos_fuera <- orden_fallo[seq_len(n_remove)]

  g_agregado_reducido <- delete_vertices(g_agregado, nodos_fuera)
  comp_red <- components(g_agregado_reducido)
  gcc_red <- if (length(comp_red$csize) == 0) 0 else max(comp_red$csize)

  tibble(
    frac_eliminada = f,
    tamano_gcc_rel = gcc_red / tamano_gcc_inicial
  )
})

resultado_robustez %>% 
  ggplot(aes(x = frac_eliminada, y = tamano_gcc_rel)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  labs(
    title = "Robustez de la red agregada ante fallos en la capa METRO",
    x = "Fracción de nodos eliminados en METRO",
    y = "Tamaño relativo de la componente gigante agregada"
  )

Una caída rápida del tamaño de la componente gigante indica alta vulnerabilidad ante fallos dirigidos en nodos centrales de la capa metro.

7 Conclusiones

En esta memoria se ha construido y analizado una red multicapas de transporte para Madrid utilizando exclusivamente los datos proporcionados:

  • Las paradas se han agrupado en zonas de intercambio mediante clustering espacial.
  • Se ha definido una red multicapa con capas para bus, metro, bicimad y parking.
  • Se ha calculado la centralidad de las zonas en cada capa y un indicador de versatilidad de grado que muestra qué zonas son importantes en varios modos a la vez.
  • Se ha generado una visualización multicapa con muxViz que hace explícita la estructura en capas.
  • Se ha estudiado la robustez estructural de la red agregada frente a fallos dirigidos en nodos centrales de la capa metro.

A partir de aquí, se podrían ampliar los análisis:

  • Incorporando información de líneas reales (no solo proximidad espacial).
  • Usando otras medidas de centralidad (betweenness, eigenvector, PageRank).
  • Implementando modelos de cascadas más cercanos a la literatura de redes interdependientes.

El código está preparado para ejecutarse directamente (asumiendo que los CSV y los paquetes necesarios están disponibles) y centrado únicamente en los datos reales de Madrid, tal y como se pedía en el enunciado.